WARNING - this work is stolen!! I have adapted this from a repository on GitHub from the wonderfully talented Milos Agathon. All credit to Milos! What a boss - really awesome stuff. Also of note is the image used for the blog. That is Cotey Gallagher… I hope she doesn’t sue me. https://www.linkedin.com/pulse/how-crazy-would-could-really-rain-cats-dogs-cotey-gallagher/
Really interested in quantifying weather data for specific areas that we are working…. Here is a first start.
First thing we will do is load our packages. If you do not have the packages installed yet you can change the update_pkgs param in the yml of this file to TRUE. Using pak is great because it allows you to update your packages when you want to.
Code
pkgs <-c("here","fs","pRecipe","giscoR","terra","tidyverse","rayshader","sf","classInt","rgl")if(params$update_pkgs){# install the pkgslapply(pkgs, pak::pkg_install,ask =FALSE)}# load the pkgsinvisible(lapply(pkgs, require,character.only =TRUE))
Define our Area of Interest
Here we diverge a bit from Milos version as we are going to load a custom area of interest. We will be connecting to our remote database using Poisson Consulting’s fwapgr::fwa_watershed_at_measure function which leverages the in database FWA_WatershedAtMeasure function from Simon Norris’ wonderful fwapg package. We use a blue line key and a downstream route measure to define our area of interest which is the Neexdzii Kwa (a.k.a Upper Bulkley River) near Houston, British Columbia.
Uniquely identifies a single flow line such that a main channel and a secondary channel with the same watershed code would have different blue line keys (the Fraser River and all side channels have different blue line keys).
A downstream route measure is:
The distance, in meters, along the route from the mouth of the route to the feature. This distance is measured from the mouth of the containing route to the downstream end of the feature.
Code
# lets build a custom watersehed just for upstream of the confluence of Neexdzii Kwa and Wetzin Kwa# blueline keyblk <-360873822# downstream route measuredrm <-166030.4aoi <- fwapgr::fwa_watershed_at_measure(blue_line_key = blk, downstream_route_measure = drm) |> sf::st_transform(4326)
# let's create our data directorydir_data <- here::here('posts', params$post_dir_name, "data")fs::dir_create(dir_data)
To actually download the data we are going to put a chunk option that allows us to just execute the code once and update it with the update_gis param in our yml. We will use usethis::use_git_ignore to add the data to our .gitignore file so that we do not commit that insano enourmouse file to our git repository.
Now we read in our freshly downloaded .nc file and clip to our area of interest and transform the data into a dataframe. We need to remove the data from 2023 because it is only for January as per the .nc filename mswep_tp_mm_global_197902_202301_025_yearly.nc.
Code
# get the name of the file with a .nc at the endnc_file <- fs::dir_ls(dir_data, pattern =".nc")# list.files()mswep_data <- terra::rast( nc_file) |>terra::crop( aoi)names(mswep_data) <-1979:2023mswep_df <- mswep_data |>as.data.frame(xy =TRUE) |> tidyr::pivot_longer(!c("x", "y"),names_to ="year",values_to ="precipitation" ) |># 2023 is not complete so we remove it dplyr::filter(year !=2023)
Plot the Precipitation Data by Year
First thing we do here is highjack the plot function from Milos.
Code
theme_for_the_win <-function(){theme_minimal() +theme(axis.line =element_blank(),axis.title.x =element_blank(),axis.title.y =element_blank(),axis.text.x =element_blank(),axis.text.y =element_blank(),legend.position ="right",legend.title =element_text(size =11, color ="grey10" ),legend.text =element_text(size =10, color ="grey10" ),panel.grid.major =element_line(color =NA ),panel.grid.minor =element_line(color =NA ),plot.background =element_rect(fill =NA, color =NA ),legend.background =element_rect(fill ="white", color =NA ),panel.border =element_rect(fill =NA, color =NA ),plot.margin =unit(c(t =0, r =0,b =0, l =0 ), "lines" ) )}breaks <- classInt::classIntervals( mswep_df$precipitation,n =5,style ="equal")$brkscolors <-hcl.colors(n =length(breaks),palette ="Temps",rev =TRUE)
Pretty cool. Interesting to see really wet and dry years in the last 20 years or so such as the wet ones in 2004, 2007, 2011 and 2020 and the dry ones in 2000, 2010, 2014, 2021 and 2022. The contours on the maps are really interesting as they show the gradients which generally run west to east - but occasionally run south to north.
Average Precipitation
Now we will average all the years together to get an average precipitation map.
Compare the Average Precipitation to a Specific Year
We often talk about a “dry” year or a “wet” year. Let’s compare the average precipitation to a specific year. We will build a little function to do this so that we can easily compare any year to the average.